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Foams have unique rheological properties that range from solid-like to fluid-like. We study two- 
^\ ' dimensional non-coarsening foams of different disorder under shear in a Monte Carlo simulation, 

0^ ■ using a driven large-Q Potts model. Simulations of periodic shear on an ordered foam show several 

0^ ' different response regimes. At small strain amplitudes, bubbles deform and recover their shapes elas- 

tically, and the macroscopic response is that of a linear elastic cellular material. For increasing strain 
r^ ■ amplitude, the energy-strain curve starts to exhibit hysteresis before any topological rearrangements 

D ' occur, indicating a macroscopic viscoelastic response. When the applied strain amplitude exceeds 

pL^ , a critical value, the yield strain, topological rearrangements (Tl events) occur, the foam starts to 

flow, and we observe macroscopic irreversibility. We find that the dynamics of topological rear- 
rangements depend sensitively on the structural disorder. Structural disorder decreases the yield 
strain; sufficiently high disorder changes the macroscopic response of a foam from a viscoelastic solid 
C^_^ ' to a viscoelastic fluid. This wide-ranging dynamical response and the associated history effects of 

O ' foams result from avalanche-like Tl events. The spatio-temporal statistics of Tl events Tl do not 

. . display long-range correlations for ordered foams or at low shear rates, consistent with experimental 

observations. As the shear rate or structural disorder increases, the topological events become more 
correlated and their power spectra change from that of white noise toward 1// noise. Intriguingly, 
the power spectra of the total stored energy also exhibit this 1// trend. 
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O PACS. 83.70.Hq, 82.70.Rr, 02.70.Lq, 46.60.Cn. 

I. INTRODUCTION 

In addition to their wide-spread industrial importance pi , foams provide significant clues to the rheology of other 

complex fluids, such as emulsions, colloids and polymer melts, because we can observe their structures directly. The 

fvj ', topological structures and the dynamics studied here also occur in other cellular materials, such as biological tissues 

f**) ■ and polycrystalline alloys. One of the most remarkable and technologically relevant features of foams is the range 

^\ ' of mechanical properties that arises from their structure. For sufficiently small stress, foams behave like a solid 

On ' and are capable of supporting static shear stress. For large stress, foams flow and deform arbitrarily like a fluid. 

■4-^ However, we do not yet fully understand the relationship between the macroscopic flow properties of foams and their 

^ , microscopic details, e.g. liquid properties, topological rearrangements of individual bubbles and structural disorder. 

^ ' Constructing a full multiscale theory of foam rheology is challenging. Foams display multiple length scales with many 

competing time scales, memory effects {e.g. the hysteresis discussed in Sec. IV below), and slow aging punctuated 

by intermittent bursts of activity {e.g. the avalanches of Tl events discussed in Sec. V below), all of which severely 

Q limit their predictability and control. These problems are intriguing both from an applied and from a fundamental 

O perspective — they provide beautiful concrete examples of multiscale materials, where structure and ordering at the 

L| microscale, accompanied by fast and slow time scales, can lead to a highly nonlinear macroscopic response. Here 

• i-H , we study the relation between the microscopic topological events and the macroscopic response in two-dimensional 

^^ ' non-coarsening foams using a driven extended large-Q Potts model. 

JH I In foams, a small volume fraction of fluid forms a continuous network separating gas bubbles ||l|. The bubble 

- - ■ shapes can vary from spherical to polyhedral, forming a complex geometrical structure insensitive to details of the 
liquid composition or the average bubble size B . Because of the complexity of describing the network of fllms and 
vertices in three-dimensional foams, most studies have been two-dimensional. In two-dimensional foams free of stress, 
all vertices are three-fold and the walls connecting them meet at 120° angles. Minimization of the total bubble wall 
length dictates that a pair of three-fold vertices is energetically more favorable than a four-fold vertex. Therefore, 
topology and dynamics are intimately related, with the dominance of three-fold vertices resulting from considerations 
of structural stability in the presence of surface tension. When shear stress is present, a pair of adjacent bubbles can be 
squeezed apart by another pair (Fig. 1), leading to a Tl switching event |^. This local but abrupt topological change 
results in bubble-complexes rearranging from one metastable configuration to another. The resulting macroscopic 
dynamics is highly nonlinear and complex, involving large local motions that depend on structures at the bubble 
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scale. The spatio-temporal statistics of Tl events is fundamental to the plastic yielding of two-dimensional liquid 
foams. 

The nonlinear and collective nature of bubble rearrangement dynamics have made analytical studies difficult, except 
under rather special assumptions. Computer simulations can therefore provide important insights into the full range 
of foam behavior. Previous studies in this field can be categorized through their use of constitutive, vertex, center or 
bubble models. 

The constitutive models have evolved from the ideas of Prud'homme Q and Princen Q. They modeled foam as a 
two-dimensional periodic array of hexagonal bubbles where Tl events occur instantaneously and simultaneously for 
the entire foam. Khan and Armstrong M further developed the model to calculate the detailed force balance at the 
films and vertices, and studied the stress-strain relationships as a function of hexagon orientation, liquid viscosity 
and liquid fraction. Reinelt and Kraynik H] extended the same model to study a polydisperse hexagonal foam and 
derived explicit relations between stress and strain tensors. While analytical calculations exist only for periodic 
structures or for linear response, foams are naturally disordered with an inherent nonlinear response. Treating the 
foam as a collection of interacting vertices, vertex models studied the eff'ect of stress on structure and the propagation 
of defects in foams with zero liquid fraction (i.e. dry foam) 1^. Okuzono and Kawasaki [g| studied the effect of 
finite shear rate by including the force on each vertex, a term which depends on the local motion and is based on 
the work of Schwartz and Princen ||9[. They predicted avalanche-like rearrangements in a slowly driven foam, with 
a power law distribution of avalanche size vs. energy release, characteristic of self-organized criticality Durian's 
|lO| , [ll| "bubble" model, treating bubbles as disks connected by elastic springs, measured foam's linear rheological 
properties as a function of polydispersity and liquid fraction. He found similar distributions for the avalanche-like 
rearrangements with a high frequency cutoff. Weaire et al. [ [l2[ , using a center model based on Voronoi construction 
from the bubble centers, applied extensional deformation and bulk shear to a two-dimensional foam. They concluded 
that avalanche-like rearrangements are possible only for wet foams, and that topological rearrangements can induce 
ordering in a disordered foam. A review by Weaire and Fortes [O] includes some computer models of the mechanical 
and rheological properties of liquid and solid foams. However, few models have attempted to relate the structural 
disorder and configuration energy to foam rheology Only recently, SoUich et al. [|l4| , studying mechanisms for storing 
and dissipating energy, emphasized the role of both structural disorder and metastability in the rheology of soft glassy 
materials, including foams. Langer and Liu |15[| , using a bubble model similar to Durian's, found that the randomness 
of foam packing has a strong effect on the linear shear response of a foam. One of the goals of our study is to quantify 
the extent of metastability by measuring hysteresis, and relate the macroscopic mechanical response to microscopic 
bubble structures. 

Experiments have measured the macroscopic mechanical properties of three-dimensional foams. But due to the 
difficulty of direct visualization in three-dimensional foams, no detailed studies of rearrangements exist. Khan et 
al. [ [16| applied bulk shear to a foam trapped between two parallel plates and measured the stress-strain response, as 
well as the yield strain as a function of liquid fraction. Princen and Kiss p7| , applying shear in a concentric cylinder 
viscometer [i.e. boundary shear), determined the yield stress and shear viscosity of highly concentrated water/oil 
emulsions. Recently, with the help of diffusing wave spectroscopy (DWS), experiments by Gopal and Durian on three- 
dimensional shaving creams showed that the rate of rearrangements is proportional to the strain rate, and that the 
rearrangements are spatially and temporally uncorrelated |^; Hohler et al. [ |l8| found that under periodic boundary 
shear, foam rearrangements cross from a linear to nonlinear regime; Hebraud et al. [ [L9[ in a similar experiment on 
concentrated emulsions, found that some bubbles follow reversible trajectories while others follow irreversible chaotic 
trajectories. However, none of these experiments has directly observed changes in bubble topology. Dennin and 
Knobler [po| performed a bulk shear experiment on a monolayer (2D) Langmuir foam and counted the number of 
bubble side-swapping events. Unfortunately, limited statistics rendered their results difficult to interpret. 

In an attempt to reconcile the different predictions of different models and experiments, we use a Monte Carlo 
model, the extended large-Q Potts model, to study foam rheology. The large-Q Potts model has successfully modeled 
foam structure, coarsening and drainage |£l|,^, capturing the physics of foams more realistically than other models. 
Here we extend the model to include the application of shear to study the mechanical response of two-dimensional 
foams under stress. 

This paper is organized as follows: Sec. II presents our large-Q Potts model; Sec. HI contains a description of 
simulation details; Sec. IV presents results on hysteresis; Sec. V discusses the dynamics and statistics of Tl events; 
Sec. VI discusses structural disorder and Sec. VII contains the conclusions. 



II. MODEL 

The great advantage of our extended large-Q Potts model is its simplicity. The model is "realistic" in that the 
position and diffusion of the walls determine the dynamics, as they do in real foams and concentrated emulsions. 
Previous models |H,|ri],^ were based on different special assumptions about the energy dissipation. Since the energy 
dissipation is poorly understood and also hard to measure in experiments, the exact ranges of validity for these models 
are not clear. Not surprisingly, these models lead to conflicting predictions, e.g. for the distribution of avalanche-like 
rearrangements (Sec. V). None of these models alone captures the full complexity of real foams. 

The extended large-Q Potts model, where bubbles have geometric properties as well as surface properties, is not 
based on any a priori energy dissipation assumption. In addition, it has the advantage of simultaneously incorporating 
many interactions, including temperature effects, for foams with arbitrary disorder and liquid content P3| . 

Both the film surface properties and the geometry of bubbles are fundamental to understanding foam flow. The 
contact angle of walls between vertices indicates whether the structure is at equilibrium, corresponding to minimizing 
the surface energy. In a real evolving pattern, the equilibrium contact angle occurs only for slow movements during 
which the vertices remain adiabatically equilibrated. Whenever a topological rearrangement (a Tl event) of the 
pattern occurs, the contact angles can be far from their equilibrium values. The walls then adjust rapidly, at a 
relaxation rate depending on the effective foam viscosity, to re-establish equilibrium. The same holds true for the 
other possible topological change, the disappearance of a bubble, a T2 event S. However, disappearance only occurs 
in foams that do not conserve bubble number and area, which we do not consider in this study. A difficulty in 
two-dimensional foams is that the effective viscosity depends primarily on the drag between the Plateau borders and 
the top/bottom surfaces of the container, not the liquid viscosity. Container chemistry, surfactant properties, and 
foam wetness all change the effective viscosity. Thus even in experiments, the effective viscosity is not equivalent to 
the liquid viscosity and not possible to derive from liquid viscosity. We define the equilibrium contact angle so that 
any infinitesimal displacement of the vertex causes a second-order variation of the surface energy, while during a Tl 
event the energy must vary macroscopically over a small but finite coherence length, typically the rigidity length of 
a bubble. In our simulations, a bubble under stress can be stretched or compressed up to 60% of its original length, 
while conserving its area. 

In a center model based on the Voronoi construction (see, e.g., |1J]), the coherence length of a bubble is comparable 
to its diameter. Contact angles are given correctly at equilibrium but approach and remain near 90° during a Tl 
event, since the centers are essentially uninfluenced by topological details such as the difference between a four-fold 
vertex and a pair of three-fold vertices. 

In a vertex model (see, e.g., Iq]), the walls connecting the vertices adiabatically follow an out-of-equilibrium, slowly- 
relaxing vertex. In such a model, the walls are constrained to be straight and vertices typically have arbitrary angles. 
In essence, the deviation of the vertex angles from the equilibrium value represents the integrated curvature of the 
bubble walls. Because of their unphysical representation of contact angles, pure vertex models with straight walls 
cannot handle Tl events correctly. 

The extended large-Q Potts model avoids these limitations: walls are free to fluctuate, which is not true in vertex 
models; and the contact angles during a Tl event are correct, which is not true in center models. A further advantage 
of our extended large-Q Potts model is that it allows direct measurement of Tl events. The other models cannot 
directly count Tl events. Instead, they quantify rearrangement events by their associated decreases in energy. We 
will discuss later this energy decrease is not always directly proportional to the number of Tl events. Our model 
therefore delivers accurate information about individual Tl events as well as the averaged macroscopic measures such 
as total bubble wall length, thereby allowing new insights into the connection between microscopic foam structure 
and macroscopic mechanical response. 

Before describing the details of the Potts model, we should first mention its major limitations. Viscosity is one 
of the basic physical properties of foams, but it is not easily specified a priori in the Potts model. Although we 
can extract the effective viscosity and the viscoelasticity of foams from simulations, we lack a clear quantitative 
description of the foam viscosity in Potts model simulations and how it relates to the effective and liquid viscosities 
of a two-dimensional foam. However, our ignorance about simulation viscosity is equivalent to our ignorance about 
experimental two-dimensional foam viscosity. Quantitative experiments will help to separate the roles of the Plateau 
borders, fluid viscosity, and topological rearrangements in determining the effective foam viscosity. A second possible 
limitation is the size effect due to lattice discretization. We show in Sec. HI that this problem does not invalidate 
our simulations. A third drawback is that the Monte Carlo algorithm results in uncertainties in the relative timing of 
events on the order of a few percent of a Monte Carlo step. While this uncertainty is insignificant for well separated 
events, it can change the measured interval between frequent events. 

The extended large-Q Potts model treats foams as spins on a lattice. Each lattice site i = {xi,yi) has an integer 
"spin" tJi chosen from {1, . . . , Q}. Domains of like spins form bubbles, while links between different spins define the 



bubble walls (films). Thus each spin merely acts as a label for a particular bubble. The surface energy resides on 
the bubble walls only. Since the present study focuses on shear-driven topological rearrangements over many loading 
cycles, we prohibit foam coarsening by applying an area constraint on individual bubbles. In practical applications, 
foam deformation and rearrangement under stress is often much faster than gas diffusion through the walls, so 
neglecting coarsening is reasonable. The Potts Hamiltonian, the total energy of the foam, includes the surface energy 
and the elastic bulk energy: 

ij n 

where J'ij is the coupling strength between neighboring spins Ci and Uj, summed over the entire lattice. The first 
term gives the total surface energy. The second term is the area constraint which prevents coarsening. The strength 
of the constraint (F) is inversely proportional to the gas compressibility; a„ is the area of the nth bubble and An its 
corresponding area under zero applied stress. We can include coarsening by setting T to zero. 
We extend the Hamiltonian to include shear: 

n' ='H + Y,l{V^.t)x^{l-5„,a,). (2) 

i 

The new term corresponds to applying shear strain (a detailed explanation follows below) to the wall between neigh- 
boring bubbles Ui and CTj, with 7 corresponding to the strain field, {xi,yi) to the coordinate of spin ai and (1,0) is 
the direction of the strain. 

The system evolves using Monte Carlo dynamics. Our algorithm differs from the standard Metropolis algorithm: 
we choose a spin at random, but only reassign it if it is at a bubble wall and then only to one of its unlike neighbors. 
The probability of accepting the trial reassignment follows the Boltzmann distribution, namely: 

r 1 An'<o , . 

^ \ exp(-A7^'/T) AH' > ' ^^' 

where A7i' is the change in Ti' due to a trial spin flip, and T is temperature. Time is measured in units of Monte Carlo 
steps (MCS), where one MCS consists of as many spin trials as there are lattice sites. This algorithm reproduces 
the same scaling as classic Monte Carlo methods in simulations of foam coarsening, but significantly reduces the 
simulation time p4|. 

The second term in Ti' biases the probability of spin reassignment in the direction of increasing Xi (if 7 < 0) or 
decreasing Xi (if 7 > 0). From dimensional analysis of 7i', 7 has units of force, but we can interpret it as the strain field 
for the following reason: In the Potts model a bubble wall segment moves at a speed proportional to the reassignment 
probability P; in this case, 

V ex ^P, (4) 

where the prefactor follows from dimensional analysis. This shear term effectively enforces a velocity v at the bubble 
walls, therefore it imposes a strain rate on the foam. The strain e(i) is then proportional to a time integral of v, 



VW)Pdt'. (5) 

If we limit the application of this term to the boundaries of the foam, we impose a boundary shear, equivalent to 
moving the boundary of the foam with no-slip between bubbles touching the boundary and the boundary, i.e., 

loG{t) yi = yinin 

-7oG(i) yi = yinax , (6) 

otherwise 

where 70 is the amplitude of the strain field and G{t) is a normalized function of time. On the other hand, 

l = |3y^G{t), (7) 

with yi between ymin and j/maxi corresponds to applying bulk shear with the strain rate varying linearly as a function 
of position in the foam. The gradient of strain rate is the shear rate, B. The corresponding experiment would be 



similar to Dcnnin and Knobler's monolayer Langmuir foam experiment |20|: a monolayer foam (2D) on the surface of 



a liquid is sheared in a concentric Couette cell, with no-slip conditions between the bubbles and the container surface. 
In all our studies we use G{t) = 1 for steady shear, and G{t) = sin(a;t) for periodic shear. Since for steady shear, the 
strain is a constant times time, or y/jPt, plotting with respect to time is equivalent to plotting with respect to strain. 
Note that our driving in the Potts model differs from that in driven spin systems, for which a large body of literature 
addresses the dynamic phase transition as a function of driving frequency and amplitude p5[ . Our driving term acts 
on the bubble walls (domain boundaries) only, while in driven spin systems e.g. the kinetic Ising model, all spins 
couple to the driving field. The resulting dynamics differ greatly. 

III. SIMULATION DETAILS 

Experimental observations show that the mechanical responses of a foam, including the yield strain, the elastic 
moduli, and the topological rearrangements, are sensitive to the liquid volume fraction [p6[ . In particular, the 
simulations of both Durian ^^ and Weaire et al. ||l^ showed a critical liquid fraction at which a foam undergoes a 
"melting transition." Although different liquid content and drainage effects can be readily incorporated in the Potts 
model ^^ , we focus on the dynamics of topological rearrangements and do not consider liquid fraction dependence of 
flow behavior, i.e. we assume the dry foam limit in this study. Also, we ignore gas diffusion across the walls, assuming 
that bubble deformation and rearrangement are much faster than coarsening. 

The definition of time (Monte Carlo steps or MCS) is not directly related to real time, but we have made choices 
to ensure that we do not under-resolve events. A shear cycle in the periodic shear case takes about 4000 MCS. In 
our simulations, a single deformed bubble recovers on a timescale of a few MCS while the relaxation of a cluster of 
deformed bubbles takes a much longer time, on the order of 10 to 100 MCS. A Tl event by definition takes one MCS 
(the short life of a four-fold vertex) , but the viscous relaxation has to average over at least the four bubbles involved 
in the Tl event, and thus lasts much longer. 

We used periodic boundary conditions in the x direction, to mitigate finite size effects. For ordered foams under 
boundary shear, we used a 400 x 100 lattice with each bubble containing 20 x 20 lattice sites; for ordered foams 
under bulk shear, we used a 256 x 256 lattice with 16 x 16 sites for each bubble. When unstressed, all the bubbles 
are hexagons, except for those truncated bubbles touching the top and bottom boundaries. In the case of disordered 
foams, we used a 256 x 256 lattice with various area distributions. We have also performed simulations using a lattice 
of size 1024 x 1024 with 64 x 64 bubbles and a lattice of size 1024 x 1024 with 16 x 16 bubbles. The results did not 
appear to differ qualitatively. A 16 x 16 bubble has a side length around 10 lattice sites, so its smallest resolvable tilt 
angle is approximately arctan(l/10) ^ 5.7°. Had lattice effects been a problem, we would have expected a significant 
difference in the simulations with bubbles of size 64, where the smallest angle is about four times smaller. But 
increasing the simulation size from 16^ to 64^ did not lead to significant changes in the quantities we measured. Thus, 
we used bubbles of size 16^ in all the simulations reported in this paper. 

Lattice anisotropy can induce artificial energy barriers in lattice simulations. All our runs use a fourth-nearest 
neighbor interaction on a square lattice, which has a lattice anisotropy of 1.03, very close to the isotropic situation 
(lattice anisotropy of 1). 

Standard quantitative measures of cellular patterns are the topological distributions and correlations, area distri- 
butions, and wall lengths — all quantities that in principle can be measured in experiments. Since the areas are 
constrained, the evolution of the area distribution is not useful. We define the topological distribution p(n) as the 
probability that a bubble has n sides; its ?7i-th moments are /i„j = "^^ p{n){n— (n))™. The area distribution p{a) and 
its second moment /i2(a) are defined in a similar fashion for the bubble areas. We use a variety of disordered foams 
with different distributions, as characterized by their p,2{n) and /i2(a). 

In practice, we generate the initial configuration by partitioning the lattice into equal-sized square domains, each 
containing 16 x 16 lattice sites. The squares alternate offsets in every other row, so the pattern resembles a brick 
wall arranged in common bond. We then run the simulation with area constraints, but without strain, at finite 
temperature for a few Monte Carlo steps, and then decrease the temperature to zero and let the pattern relax. The 
minimization of total surface energy (and hence the total bubble wall length) results in a hexagonal pattern, the initial 
configuration for the ordered foam. For disordered initial configurations, we continue to evolve the hexagonal pattern 
without area constraints at finite temperature so that the bubbles coarsen. We monitor /i2 {n) of the evolving pattern, 
and stop the evolution at any desired distribution or degree of structural disorder. Then we relax the patterns at 
zero temperature with area constraints to guarantee that they have equilibrated, i.e. without added external strain 
or stress the bubbles would not deform or rearrange. 

For all our simulations, P = 1 (which is sufhciently large to enforce air incompressibility in bubbles) and J7ij = 3 
(except when we vary the coupling strength to change the effective viscosity of the foam) . Most of the simulations 
shown in this paper are run at zero temperature except when we study temperature effects on hysteresis, because the 



data are less noisy and easier to interpret. A finite but low temperature speeds the simulations, but does not appear 
to change the results qualitatively. 

The number of sides of a bubble is defined by its number of different neighbors. During each simulation, we keep 
a list of neighbors for each bubble. A change in the neighbor list indicates a topological change which, since bubbles 
do not disappear, has to be a Tl event. 

IV. HYSTERESIS 

We can view foam flow as a collective rearrangement of bubbles from one metastable configuration to another. We 
investigate the configurational metastability by studying hysteresis of the macroscopic response. 

Hysteresis is the phenomenon in which the macroscopic state of a system does not reversibly follow changes in an 
external parameter, resulting in a memory effect. Hysteresis commonly appears in systems with many metastable 
states due to (but not limited to) interfacial phenomena or domain dynamics. The classic example of the former is 
that the contact angle between a liquid and a solid surface depends on whether the front is advancing or retreating. 
The classic example of the latter is ferromagnetic hysteresis, in which the magnetization lags behind the change 
in applied magnetic field. In cellular materials, including foams, hysteresis can have multiple microscopic origins, 
including stick-slip interfacial and vertex motion, local symmetry-breaking bubble rearrangement (Tl events), and 
the nucleation of new and annihilation of old cells. In all of these, noise and disorder play an intrinsic role in selecting 
among the many possible metastable states arising when the foam is driven away from equilibrium. By focusing on 
non-coarsening foams, we rule out nucleation and annihilation as sources for hysteresis. Our foam is therefore an ideal 
testing ground for improving our understanding of hysteresis as it arises from local rearrangements and interfacial 
dynamics. 

In accordance with |g[ , we define the quantity: 

0^^0(1-<5..,,J, (8) 

as the total stored elastic energy. Here sites i,j are neighbors, summation is over the whole lattice, and 9 is the wall 
thickness that we choose to be 1 in all our simulations (dry foam limit) . Thus gives essentially the total bubble wall 
length, which differs by a constant, namely the surface tension, from the total surface energy. In zero temperature 
simulations, the area constraint is almost always satisfied so that small fluctuations in areas contribute only 10~^ of 
the total energy. Thus we can neglect the elastic bulk energy of the bubbles, and assume that the total foam energy 
resides on the bubble walls only, i.e. all forces concentrate at the bubble walls. We can calculate values of the averaged 
stress by taking numerical derivatives of the total surface energy with respect to strain [ |l2| . However, the calculation 
via derivatives is not suitable for foams undergoing many topological changes, since the stored elastic energy changes 
discontinuously when topological rearrangements occur. The alternative is to calculate stress directly, as given in 
pTj , by the sum of forces acting on the bubble walls, which locally is proportional to the wall length change of a 
bubble. Because forces on the bubble walls in Potts model foams are not well characterized, we limit our discussions to 
energy-strain relationships. The more rigorous definition of strain involves the definition of a mesoscopic lengthscale 
corresponding to a cluster of bubbles, over which the effects of bubble wall orientation and bubble deformation can 
be averaged. In [gj, the average stress tensor, defined as cr = l/^X](i j) \'''ij\^ij^i]^ with A the total area of the foam 
and Tij the distance between two neighboring vertices, is directly related to via — Tr{a). Hereafter, we present 
our (j) data as 4>{t)/(p{0) to scale out differences due to initial configurations. 

A. Hysteresis in ordered foams 

The simplest perturbation which induces topological rearrangements is boundary shear on an ordered foam. In 
this case we can confine the deformation to the bubbles touching the moving boundaries, and easily locate all the 
Tl events. As the applied boundary shear increases, the bubbles touching the boundaries distort, giving rise to a 
stored elastic energy. We show snapshots of the pattern in Fig. 2(a). When a pair of vertices come together to 
form a four-fold vertex, the numbers of sides changes for the cluster of bubbles involved. Different shades of gray 
in Fig. 2(a) reflect the topologies of the bubbles. Note that a five-sided (dark grey) and a seven-sided (light grey) 
bubble always appear in pairs except during the short lifetime of a four-fold vertex (when the number of sides are 
ambiguous because of the discrete lattice). Once the strain exceeds a critical value, the yield strain, all the bubbles 
touching the moving boundaries undergo almost simultaneous rearrangements, thereby releasing stress. The stored 



elastic energy, (j), increases with time when the bubbles deform, then decreases rapidly when the bubbles rearrange. 
Stress accumulates only in the two boundary layers of bubbles, and never propagates into the interior of the foam. 
The whole process repeats periodically, due to the periodic bubble structure, as shown in Fig. 2(b), the energy-strain 
plot (as mentioned at the end of Sec. II, for steady shear, plotting time is equivalent to plotting strain). This result 
corresponds to the mechanical response obtained in the model of Khan et al. with periodic hexagonal bubbles oriented 
at zero degrees with respect to applied strain Q . 

When applying periodic shear 7(i) — 7osin(Li;i), we keep the period 2tt/u! fixed and vary the amplitude, 70. Under 
sinusoidal periodic shear, we observe three types of behavior. When the strain amplitude is small, bubbles deform and 
recover their shapes elastically when stress is released. No topological rearrangement occurs and the energy-strain 
plot is linear, corresponding to an elastic response [2^. This result agrees perfectly with the experimental result of 
DWS in ||l8|. As the strain amplitude increases, the energy-strain curve begins to exhibit a small butterfly-shaped 
hysteresis loop before any topological rearrangements occur, indicating a macroscopic viscoelastic response. If we 
keep increasing the strain amplitude, the hysteresis loop increases in size. When the applied strain amplitude exceeds 
a critical value, Tl events start occurring, and the foam starts to flow, which leads to a further change in the shape 
of the hysteresis loop. Even larger strain amplitude introduces more Tl events per period, and adds small loops to 
the "wings" of the hysteresis loop. Figure 3(a) shows the smooth transition between the three types of hysteresis in 
the energy-strain curve. 

We can adjust the viscosity of the bubble walls by changing the coupling strength J^ij . Smaller coupling strength 
corresponds to lower viscosity. Similar transitions from elastic, to viscoelastic to fluid-like flow behavior occur for 
progressively lower values of coupling strength, shown in Fig. 3(b). The phase diagram in Fig. 3(c) summarizes 
44 different simulations and shows the elastic, viscoelastic and fluid-like behavior (as derived from the hysteretic 
response) as a function of the coupling strengths jTij {i.e. viscosity) and strain amplitudes 79. A striking feature is 
that the boundaries between these regimes appear to be linear. Figure 3(d) shows the effect of finite temperature on 
the energy-strain curves. With progressively increasing temperature, noise becomes more dominant and eventually 
destroys the hysteresis loop. This result implies diminished metastability at finite temperature. However it does not 
seem to change the trend in mechanical response. 

A more conventional experiment is the application of bulk shear |p|, pl]jl7|j2C| ] , with the shear strain varying linearly 
as a function of the vertical coordinate, from 70 at the top of the foam to —70 at the bottom. In our bulk shear 
simulations with an ordered foam, the energy-strain relationship has two distinct behaviors depending on the shear 
rate. At small shear rates, a "sliding plane" develops in the middle of the foam. As shown in Fig. 4(a), non- hexagonal 
bubbles appear only at the center plane. The energy-strain curve, shown in Fig. 4(b), therefore, resembles that 
for boundary shear on an ordered foam. The energy curve in Fig. 4(b) also shows that the baseline of energy is 
larger and that the decrease in amplitude of the energy due to Tls is smaller than in Fig. 2(b), because bulk shear 
induces a more homogeneous distribution of distortion and thus of stored elastic energy. Again the periodic structure 
of the bubbles cause the periodicity of the curve, reminiscent of the shear planes observed in metallic glasses in the 
inhomogeneous fiow regime, where stress-induced rearrangement causes plastic deformation p9| . 

At high shear rates, the ensemble of Tl events no longer localize in space [Fig. 5(a)]. Non-hexagonal bubbles 
appear throughout the foam. The energy-strain curve, shown in Fig. 5(b), is not periodic but rather smooth; beyond 
the yield point, the bubbles constantly move without settling into a metastable configuration and correspondingly, 
the foam displays dynamically induced topological disorder. The transition between these two regimes, localized and 
nonlocalizcd Tl events, occurs when the shear rate is in the range 1 x 10"^ < |/3| < 5 x 10~^. This transition can be 
understood if we look at the relaxation time scale of the foam. Due to surface viscous drag and geometric confinement 
of other bubbles, the relaxation time for a deformed bubble in a foam is on the order of 10 MCS. For a shear rate 
/3 = 5 X 10~^, P~^ is of the same order as the relaxation time. Thus for shear rates above the natural internal 
relaxation time scale, the macroscopic response changes from jagged and piece-wise elastic to smooth and viscous 
response, as observed in fingering experiments in foams [ pO| . 

B. Hysteresis in disordered foams 

In a disordered foam, bubbles touching the moving foam boundary have different sizes. Boundary strain causes 
different bubbles to undergo Tl rearrangements at different times. Stress no longer localizes in (sliding) boundary 
layers, but propagates into the interior [Fig. 6(a)]. The yield strain is much smaller. When the size distribution of the 
foam is broad, the linear elastic regime disappears, since even a small strain may lead to topological rearrangements 
of small bubbles. In other words, with increasing degree of disorder, the yield strain decreases to zero and the foam 
changes from a viscoelastic solid to a viscoelastic fluid. We show an example of such viscoelastic fluid behavior in Fig. 
6(b) for a random foam, which shows no energy accumulation, namely, its yield strain is zero. The foam deforms and 



yields like a fluid upon application of the smallest strain. 

Under a periodic shear, the stored energy increases during an initial transient period but reaches a steady state 
after a few periods of loading. Energy-strain plots show hysteresis due to topological rearrangements similar to those 
in ordered foams, but as the degree of disorder increases, the corresponding elastic regime shrinks and eventually 
disappears. 

Rearrangement events in a disordered foam under bulk shear at a low shear rate [snapshots shown in Fig. 7(a)] 
correspond to those in an ordered foam at a high shear rate. The rearrangements are discrete and avalanche-like, 
resembling a stick-slip process, or adding sand slowly to a sandpile. However, at sufficiently high shear rate all the 
avalanches overlap and the deformation and rearrangements are more homogeneous and continuous, as in a simple 
viscous liquid. Figure 7(b) shows the typical energy-strain curve of a disordered foam under steady bulk shear. 

Note that in all our hysteresis plots, the energy-strain curves cross at zero strain, indicating no residual stored 
energy at zero strain. This crossing is an artifact of our definition of energy, which ignores angular measures of 
distortion, i.e. the total bubble wall length does not distinguish among the directions in which the bubbles tilt. A 
choice of stress definition which included angular information would show some residual stress at zero strain, but 
would not affect the results reported here. 

V. Tl AVALANCHES 

In both experiments pi^,p3l and our simulations, the contact angles of the vertices remain close to 120° until two 
vertices meet. The applied strain rate determines the rate at which vertices meet. The resulting four-fold vertex 
rapidly splits into a vertex pair, recovering 120° contact angles, at a rate determined by the viscosity. This temporal 
asymmetry in the Tl event contributes to the hysteresis. 

In vertex model simulations, sudden releases of energy occur once the applied shear exceeds the yield strain 1^. The 
event size n, i.e. the energy release per event in the dry foam limit, follows a power-law distribution: p{n) ^ n^^, 
Durian ||ll[ found a similar power-law distribution in his bubble model, with an additional exponential cutoff for large 
events. Simulations of Weaire et al. [|l^,^, however, suggested that power-law behavior only appeared in the wet 
foam limit. Experiments, on the other hand, have never found system-wide events or long range correlations among 
events p,po|. One of our goals is to reconcile these different predictions. 

These differences may result from the use of energy release as proxy for topological changes, rather than enumerating 
actual events, as well as the assumption of a linear relation between jumps in the stored elastic energy and the number 
of Tl events, namely, d(j)/dt = cN, where N is the number of Tl events and c is a constant. A drastic drop in the total 
bubble wall length indicates a large number of Tl events. However, in a disordered foam, all Tl events are not equal, 
since they do not all release the same amount of stored elastic energy. The energy released during a Tl event scales 
as the bubble perimeter, i.e. smaller bubbles release less energy. Hence smaller bubbles undergo more Tl events. 
Moreover, a Tl event is not strictly local, but deforms its neighborhood over a certain finite range, as demonstrated 
by Tl manipulations in magnetic fluid foam experiments Q|. Therefore, the number of Tl events is not always 
directly proportional to the decrease in total bubble wall length. Thus we cannot compare the energy dissipation 
and Tl events directly. Furthermore, the mechanisms of energy dissipation differ in these models. Kawasaki et al. M 
included the dissipation due to the flow of liquid out of the Plateau borders; Durian |l^ considered only the viscous 
drag of the liquid, while Weaire et al. |l3| modeled an equilibrium calculation involving quasistatic steps in the strain 
that do not involve any dissipation. In our model, the evolution minimizes the total free energy naturally. To avoid 
ambiguities, we directly count Tl events in addition to tracking energy. 

The avalanche-like nature of rearrangements appears in the sudden decreases of the total elastic energy as a function 
of time. Figure 4(b) shows the relation between energy and the number of Tl events in an ordered foam under steady 
bulk shear for a small strain. The stored energy increases almost linearly until the yield strain is reached. The 
avalanches are well separated. Every cluster of Tl events corresponds to a drastic decrease in the stress, and the 
periodicity is due to the ordered structure of the foam. At a higher shear rate [Fig. 5(b)], the yield strain remains 
almost the same, but the avalanches start to overlap and the energy curve becomes smoother. In the sandpile analogy, 
instead of adding sand grains one at a time and waiting until one avalanche is over before dropping another grain, 
the grains accumulate at a constant rate and the avalanches, large and small, overlap one another. A sufficiently 
disordered foam may not have a yield strain [Fig. 6(b)]; Tl events occur at the smallest strain. The foam flows as a 
fluid without going through an intermediate elastic regime. 

To study the correlation between Tl events, we consider the power spectrum of N(t), the number of Tl events at 
each time step, 

PNif) ^ fdt f dTe-'f^N{t)N{t + r), (9) 



where / is the frequency with unit MCS~^ . Figure 8(a) shows typical power spectra of the time series of Tl events in 
an ordered foam under bulk shear. At a shear rate /3 = 0.01, the Tl events show no power law. The peak at ~ 10"'^ is 
due to the periodicity of bubble structure in an ordered foam when a "sliding plane" develops. At shear rate [3 = 0.02, 
the spectrum resembles that of white noise. As the shear rate increases to /3 = 0.05, the power spectrum develops 
a power law tail at the low frequency end, with an exponent very close to 1. In a disordered foam, with increasing 
shear rate, the spectra for the Tl events gradually change from completely uncorrelated white noise to 1// at higher 
shear rates. By 1//, we mean any noise of power spectrum S{f) ^ /~" where < a < 2 or near 1, i.e. intermediate 
between Brownian noise (a = 2) and white noise (a = 0). 

These power spectra suggest that the experimental results for Tl events as reported in [^ 20| correspond to a low 



shear rate, with no long-range correlation among Tl events. Structural disorder introduces correlations among the 
events. Power-law avalanches do not occur in ordered hexagonal cells at low shear rate, where rearrangements occur 
simultaneously. At a high shear rate, when the value of 7~^ is comparable to the duration of rearrangement events, the 
bubbles move constantly. The foam behaves viscously, since rearrangements are continuously induced before bubbles 
can relax into metastable configurations which can support stress elastically. At these rates, even an initially ordered 
structure behaves like a disordered one, as shear destroys its symmetry and periodicity. 

In a disordered foam, whenever one Tl event happens, the deformed bubbles release energy by viscous dissipation 
and also transfer stress to their neighboring bubbles, which in turn are more likely to undergo a Tl switch. Thus, 
Tl events become more correlated. Shown in Fig. 8(b), the power spectra change from that of white noise toward 
1// noise. When the first sufficiently large region to accumulate stress undergoes Tl events, it releases stress and 
pushes most of the rest of the bubbles over the brink, causing an "infinite avalanche": some bubbles switch neighbors, 
triggering their neighbors to rearrange (and so on), until a finite fraction of the foam has changed configuration, 
causing a decrease in the total stored energy, mimicking the cooperative dynamic events in a random field Ising model 
p4| . We never observe system- wide avalanches as claimed in the vertex model simulations S , agreeing with Durian's 
simulations pl| and Dennin et al.'s experiments [ po| . For even greater disorder, the bubbles essentially rearrange 
independently, provided spatial correlations for area and topology are weak. Pairs of bubbles switch as the strain 
exceeds their local yield points. Although more frequent, the avalanches are small, without long correlation lengths. 
Figure 8(c) shows the power spectra for Tl events for a highly disordered structure with fJ,2{n) = 1.65. We observe no 
power law behavior, even at high shear rates. Thus a highly disordered foam resembles a homogeneous but nonlinear 
viscous fluid. 

Over a range of structural disorder the topological rearrangement events are strongly correlated. The question 
naturally arises whether the transition between these correlated and uncorrelated regimes is sharp or smooth, and 
what determines the transition points. We are currently carrying out detailed simulations involving different structural 
disorder to study this transition. 

Previous simulations |ll[| a nd experiments |20[| measured N, the average number of Tl events per bubble per unit 
shear, and concluded that TV was independent of the shear rate. Our simulation results in three different foams with 
shear rates covering two orders of magnitude, however, disagree. As shown in Fig. 9, our data indicate that N 
depends sensitively on both the polydispersity of the foam and the shear rate. Only at large shear rates does N seem 
to be independent of the shear rate, which might correspond to the above-mentioned experiments. 

The avalanches and 1// power spectra resemble a number of systems with many degrees of freedom and dissipative 
dynamics which organize into marginally stable states Pq] . Simple examples include stick-slip models, driven chains 
of nonlinear oscillators, and sandpilc models. In sandpile models, both the energy dissipation rate (total number 
of transport events at each time step) and the output current (the number of sand grains leaving the pile) show 
power-law scaling in their distributions. In particular, if the avalanches do not overlap, then the power spectrum of 
the output current follows a power law with a finite size cutoff pq |. The l//-type power spectra result from random 
superposition of individual avalanches |37[ . 

If the analogy with sandpiles holds, we should expect the power spectra of the time derivative difi/dt, of the stored 
energy, i.e. the energy change at every time step, to be l//-like, and thus the power spectra of (j) to be /~^. However, 
in our simulations d(j)/ dt does not show l//-type broadband noise. Figures 10 (a-c) show the corresponding power 
spectra for (j) from the same simulations as Fig. 8, which are obviously not /~^, i.e. the topological rearrangements 
are not in the same universality class as sandpiles. In particular. Fig. 8(c) shows a complicated trend: the power 
spectrum changes from a small slope at shear rate (3 = 0.001 to /^"'^ spanning over 4 decades at /3 = 0.005. But as the 
shear rate increases, the power law disappears. Instead, a flat tail develops at high frequencies due to Gaussian noise. 
Other different slopes appear over different regimes of different sizes, indicating the existence of multiple time-scales 
and length-scales. We will further explore the implications of these spectra for elsewhere [p8[ . 



VI. EFFECTS OF STRUCTURAL DISORDER 

As structural disorder plays an important role in mechanical response, we study the effect of disorder on the yield 
strain and the evolution of disorder in foams under shear. We define the yield strain, at which the first Tl avalanches 
occur, as the displacement at the top boundary of the foam divided by half the height of the foam (since the zero 
strain is in the middle of the foam) rescaled by the average bubble width. 

Figure 11(a) shows the yield strain as a function of shear rate (3 for different foam disorders. We find that for an 
ordered foam at low shear rates, when a sliding plane occurs in the middle of the foam, the yield strain is independent 
of shear rate. We expect this independence because Tl events occur almost simultaneously in the sliding plane, and 
the bubble size determines the yield strain. At high shear rates, Tl events distribute more homogeneously throughout 
the foam which lowers the yield strain. The upper limit for the yield strain in an ordered foam is 2/-\/3, when all 
the vertices in a hexagonal bubble array simultaneously become four-fold under shear. The nucleation of topological 
defects (5-and- 7-sided bubble pairs) and their propagation in foams lower the yield strain. But the yield strain does 
not reach zero even at a very high shear rate of /3 = 0.05. An ordered foam remains a solid with finite yield strain. 
For a disordered foam, the yield strain is lower for higher shear rates; and at the same shear rate, the yield strain 
decreases drastically to zero as disorder increases — the foam changes from a viscoelastic solid to a viscoelastic fluid. 

The most commonly used measure for topological disorder is the second moment of the topological distribution, 
/i2(?T.). During diffusional foam coarsening, the topological distribution tends to a stationary scaling form a,nd u^in 



assumes a roughly constant value. Experiments on soap foams with up to 10000 bubbles in the initial state [^1 32 
and early smaller simulations p9| gave a value oi U2{n) = 1.4 in the scaling regime. Other simulations showed a 
slightly lower value of yLt2 = 1.2 [}40|. Weaire et al. llj] reported shear-induced ordering, i.e. reduction oi ii2{n) with 
shearing. However, in our simulations, foams with initial ^2{n) ranging from 0.81 to 2.02 show no shear-induced 
ordering. Instead, /i2("-) increases and never decreases back to its initial unstrained value. Figure 11(b) shows the 
evolution of iJi2{n) for a variety of initial topological distributions. The difference between the simulations of Weaire 
et al. 112] and ours is not surprising. Weaire et al. applied step strain and observed the resulting equilibrated pattern. 
In our simulations, bubbles are constantly under shear, i.e. the foam is not in equilibrium. The topological disorder, 
as measured by iJ,2{n), therefore increases as the energy accumulates and decreases as the energy releases, as does the 
number of topological events, and does not necessarily settle to an equilibrium value. 

In an ordered foam at shear rate /3 = 0.01 [Fig. 4] with separated Tl avalanches, ^2{n) fluctuates in synchrony 
with the total energy, shown in Fig. 11(c). When the avalanches overlap, H2{n) fluctuates more smoothly, but almost 
always has a positive correlation with the total stored energy. 

Notice that in the energy-strain plots [Fig. 2(b), Fig. 4(b)], stored energy slowly increases over long times, because 
we continuously apply shear and the foam is always out of equilibrium. Bubbles do not fully recover their original 
shapes. This deformation slowly accumulates at long times. In disordered foams, topological rearrangement may 
enhance the spatial correlation of bubbles, i.e. small bubbles cluster over time, as predicted by Langer and Liu's 
bubble model [|l^. We will report results on spatial correlations for bubble topology n, and area a elsewhere p^ . 

VII. CONCLUSIONS 

We have included a driving term in the large-Q Potts model to apply shear to foams of different disorder. When 
the driving rate is too fast for the foam to relax, the system falls out of equilibrium. The mechanical response then 
lags behind the driving shear, resulting in hysteresis. Our model differs from most well studied driven spin models: 
our spins do not couple to an external field the way Ising spins couple to an oscillating magnetic field, and all action 
occurs only at the domain boundaries. 

Because of the difficulty in characterizing local stress and strain in Potts model foams, we have chosen to use 
the rescaled total bubble wall length, (p, as the order parameter for hysteresis. While the hysteresis loops reflect 
the nonlinearity and metastability of bubble configurations, it is still an open question whether we can find more 
appropriate order parameter(s) that will provide more insight into the dynamics of Tl events. Another consequence 
of this difhculty is the lack of a clear quantitative description of the viscosity in our simulated foams in terms of the 
model parameters, a difficulty mirrored in the lack of understanding of effective foam viscosity in experiments. As 
mentioned above, a fundamental problem is the lack of experimental data on the viscosity of two-dimensional foams. 
We hope that these simulations will motivate new experiments in this direction. 

The local cellular patterns characteristic of Tl events in foams are strikingly similar to the low temperature defects 
and the hexatic-square Voronoi patterns observed in two-dimensional (particle) systems, e.g. two-dimensional liquid 
crystals and colloidal suspensions, where studies have focused on the melting phase transitions [^,Q. This similarity 
lead us to try the defect description used in melting studies, namely the nearest-neighbor-bond-orientation order 
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parameter ^„exp(6i0„), where 6'„ is the angle between two neighboring bonds. However, we found it insensitive to 
the orientation change of bubble walls during Tl events and thus not a useful order parameter. How the nucleation 
and propagation of topological defects in a sheared foam relate to the nucleation and role of topological defects in the 
two-dimensional melting studies remains an interesting question. 

We have demonstrated three different hysteresis regimes in an ordered foam under oscillating shear. At small strain 
amplitudes, bubbles deform and recover their shapes elastically after stress release. The macroscopic response is that 
of a linear elastic solid. For larger strain, the energy-strain curve starts to exhibit hysteresis before any topological 
rearrangements occur, indicating a macroscopic viscoelastic response. Increasing the strain amplitude increases the 
area of the hysteresis loop. When the applied strain amplitude exceeds a critical value, the yield strain, Tl events 
occur and the foam starts to flow, and we observe macroscopic irreversibility. 

We are currently testing this observation in an experiment similar to p^ , applying periodic boundary shear to a 
homogeneous foam, and measuring the total bubble wall length directly (instead of using diffusion wave spectroscopy) 
to obtain (f>. We can directly compare this data with the predicted three distinct behaviors. The viscoelasticity of 
foams is better characterized using the complex modulus G{(jj) |14[ , which we plan to use in future investigations. 

The comparison between the mechanical responses of ordered and disordered foams provides some insight into the 
relation between local structure and macroscopic response. An ordered foam has a finite yield strain. Structural 
disorder decreases the yield strain; sufficiently high disorder changes the macroscopic response of a foam from a 
viscoelastic solid to a viscoelastic fluid. A random foam with broad topology and area distributions lacks the linear 
elastic and viscoelastic solid regimes. Any finite stress can lead to topological rearrangements of small bubbles and thus 
to plastic yielding of the foam. More detailed simulations and experiments are needed to determine the dependence 
of the yield strain on the area and topological distributions of the foam, and on the shear rates. High shear rates 
effectively introduce more topological defects into the foam, as manifested in ordered foams driven at high shear rates. 
Local topological rearrangements (the appearance of non-hexagonal bubbles) occur throughout the foam, resulting in 
more homogeneous flow behavior, as in disordered foams. 

Our simulations show that N, the average number of Tl events per bubble per unit shear, is sensitive to the area 
distribution of the foam and the shear rate. Only for a small range of shear rates do foams having similar distributions 
show similar values of TV, which may explain previous studies |0,EQ]. Our results emphasize the importance of both 
structural disorder and configurational metastability to the behavior of soft cellular materials. 

In disordered foams, the number of Tl events is not directly proportional to the elastic energy release, because each 
Tl event is non-local and different Tl events can release different amounts of energy. Therefore, we count Tl events 
and energy release separately. 

Avalanche-like topological rearrangements play a key role in foam rheology. Our simulations show that Tl events 
do not have finite long-range correlations for ordered structures or at low shear rates, consistent with experimental 
observations. As the shear rate or structural disorder increases, the topological events become more correlated. Over 
a range of disorders, the power spectra are 1//. As Hwa and Kardar pointed out, 1// noise may arise from a random 
superposition of avalanches |0. These 1// spectra suggest that avalanches of different sizes, although they overlap, 
are independent of each other. Both greater structural disorder and higher shear rate introduce a flat tail at the high 
frequency end, a signature of Gaussian noise, but do not change the exponent in the power law region. 

However, unlike the sandpile model, the power spectra of the total energy, rather than of the energy dissipation, 
show a similar trend toward 1//. One major difference between Tl avalanches and sand avalanches is that each sand 
grain carries the same energy, while each Tl event can have a different energy. A better analogy may be a "disordered 
sandpile" model, where the sand grains have different sizes or densities, and avalanches overlap. 

Avalanches of Tl events decrease the stored elastic energy, leading to foam flow. How do single Tl events contribute 
to the global response? Magnetic fluid foam experiments offer a unique opportunity to locally drive a vertex and 
force a single Tl event (or a Tl avalanche) by a well-controlled local magnetic field. We are investigating the effects 
of single Tl events using magnetic fluid foam experiments and the corresponding Potts model simulations 143]. 
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FIGURE CAPTIONS 

Figure 1. Schematic diagram of a Tl event, where bubbles a,b,c and d swap neighbors. Notice that as the edge 
between the pair of vertices shrinks, the contact angles not in contact with this edge remain 120°. 

Figure 2. An ordered foam under boundary shear: (a) Snapshots; different shades of grey encode bubble topologies 
(lattice size 256 x 256). (b) Energy-strain curve and the number of Tls presented in 50 MCS bins. 

Figure 3. Energy-strain curves for ordered foam under periodic boundary shear: (a) Numbers above the figures are 
7o. Progressively increasing shear amplitude at J^ij — 3 leads to a transition between three types of hysteresis: 70 = 1.0 
corresponds to an elastic response, 70 = 3.5 shows viscoelastic response (before any Tl event occurs), 70 = 7.0 is a 
typical response when only one Tl event occurs during one cycle of strain loading. The intermediate steps show that 
the transition between these three types is smooth, (b) Numbers above the figures are J7ij. Progressively decreasing 
liquid viscosity (increasing J7ij at 70 = 7) shows a similar transition between elastic (J7ij — 10) and viscoelastic 
{>Jij — 5) regimes, and flow due to Tl events {J'ij — 3 for one Tl event and J'ij = 1 for three Tl events during one 
strain cycle, (c) Phase diagram of hysteresis in the parameter space 70 vs. Jij ; (d) Effect of progressively increasing 
temperature T (^Jij = 3, 70 = 4). All data shown here are averaged over 10 periods. 

Figure 4. An ordered foam under bulk shear with shear rate /3 — 0.01: (a) Snapshots; shades of grey encode bubble 
topologies as in Fig. 2. (b) Energy-strain curve and the number of Tls. The magnified view in the box shows the 
correlation between stress releases and overlapping avalanches of Tl events. 

Figure 5. An ordered foam under bulk shear with shear rate /3 — 0.05: (a) Snapshots; (b) Energy-strain curve 
and the number of Tls. The magnified view in the box shows the correlation between stress releases and overlapping 
avalanches of Tl events. 

Figure 6. A disordered foam under boundary shear: (a) Snapshots; shades of grey encode bubble topologies (lattice 
size 256 x 256). (b) Energy-strain curve and number of Tl events presented in 100 MCS bins. 

Figure 7. A disordered foam under bulk shear at shear rate /3 = 0.01: (a) Snapshots; shades of grey encode bubble 
topologies (lattice size 256 x 256): (b) Energy-strain curve and number of Tl events. 

Figure 8. Power spectra of the number of Tl events: (a) An ordered foam for three shear rates /3 = 0.01, 0.02 and 
0.05 respectively; (b) A disordered foam \[i.i(n) = 0.81, /U2(a) = 7.25] for five shear rates from 0.001 to 0.05; (c) A 
very disordered foam \[i2{n) = 1.65, /i2(a) ~ 21.33] for three shear rates. 

Figure 9. Number of Tl events per unit shear per bubble as a function of shear rate for four foams: squares 
correspond to a foam of 180 bubbles, with iJ,2{n) = 1.65, /i2(a) — 21.33; stars correspond to a foam of 246 bubbles 
with iJ,2(n) ~ 1.72, /i2(a) = 15.1; triangles correspond to a foam of 377 bubbles, with ii2{n) — 1.07, /i2(a) = 2.50; and 
circles correspond to a foam of 380 bubbles, with /i2(«) = 0.95, /i2(a) ~ 2.35. The inset shows on a log-log scale that 
N varies by several orders of magnitude. 

Figure 10. Power spectra of the energy: (a) An ordered foam for three shear rates, (b) A disordered foam 
[/i2(n) = 0.81, /i2(a) = 7.25] for five shear rates, (c) A very disordered foam [ii2{N) — 1.65, /i2(a) = 21.33] for five 
shear rates. 

Figure 11. (a) Yield strain as a function of shear rate, (b) Evolution of ^2 {n) under constant bulk shear, legend 
denotes the initial /i2 (n) . (c) Evolution of /X2 (n) under steady bulk shear for an ordered foam, showing the correlation 
between the stress decreases and ii2{n); 
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